#include <stdlib.h>
#include <stddef.h>
#include "mt19937.h"


#define N 624
#define M 397
#define MATRIX_A 0x9908b0df /* constant vector a */
#define UPPER_MASK 0x80000000 /* most significant w-r bits */
#define LOWER_MASK 0x7fffffff /* least significant r bits */

#define NN 312
#define MM 156
#define MATRIX_A_64 0xB5026F5AA96619E9ULL
#define UM 0xFFFFFFFF80000000ULL /* Most significant 33 bits */
#define LM 0x7FFFFFFFULL /* Least significant 31 bits */

struct __mt19937_engine{
    unsigned int mt[N];
    unsigned mti;
};

struct __mt19937_64_engine{
    unsigned long long mt[NN];
    unsigned mti;
};

/* initializes mt[N] with a seed */
mt19937_engine *create_mt19937(unsigned s)
{
    unsigned mti;
    mt19937_engine *engine=(mt19937_engine*)malloc(sizeof(mt19937_engine));
    engine->mt[0]= s & 0xffffffff;
    for (mti=1; mti<N; mti++) {
        engine->mt[mti] =
                (1812433253 * (engine->mt[mti-1] ^ (engine->mt[mti-1] >> 30)) + mti);
        /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
        /* In the previous versions, MSBs of the seed affect   */
        /* only MSBs of the array mt[].                        */
        /* 2002/01/09 modified by Makoto Matsumoto             */
        engine->mt[mti] &= 0xffffffff;
        /* for >32 bit machines */
    }
    engine->mti=mti;
    return engine;
}

/* initialize by an array with array-length */
/* init_key is the array for initializing keys */
/* key_length is its length */
mt19937_engine *create_mt19937_by_array(const unsigned *init_key, size_t key_length)
{
    size_t i=1, j=0, k= (N>key_length ? N : key_length);
    mt19937_engine *engine=create_mt19937(19650218);
    for ( ; k; k--){
        engine->mt[i] = (  engine->mt[i] ^ ( ( engine->mt[i-1] ^ (engine->mt[i-1] >> 30) ) * 1664525 )  ) + init_key[j] + j; /* non linear */
        engine->mt[i] &= 0xffffffff; /* for WORDSIZE > 32 machines */
        i++; j++;
        if (i>=N) { engine->mt[0] = engine->mt[N-1]; i=1; }
        if (j>=key_length) j=0;
    }

    for (k=N-1; k; k--) {
        engine->mt[i] = (engine->mt[i] ^ ((engine->mt[i-1] ^ (engine->mt[i-1] >> 30)) * 1566083941))
                - i; /* non linear */
        engine->mt[i] &= 0xffffffff; /* for WORDSIZE > 32 machines */
        i++;
        if (i>=N) { engine->mt[0] = engine->mt[N-1]; i=1; }
    }
    engine->mt[0] = 0x80000000; /* MSB is 1; assuring non-zero initial array */
    return engine;
}

void destroy_mt19937(mt19937_engine *engine){
    free(engine);
}

/* generates a random number on [0,0xffffffff]-interval */
static unsigned mt19937_randuint32_nolock(mt19937_engine *engine)
{
    unsigned y;
    static const unsigned mag01[2]={0, MATRIX_A};
    /* mag01[x] = x * MATRIX_A  for x=0,1 */

    if (engine->mti >= N) { /* generate N words at one time */
        size_t kk;
        for (kk=0;kk<N-M;kk++) {
            y = (engine->mt[kk]&UPPER_MASK)|(engine->mt[kk+1]&LOWER_MASK);
            engine->mt[kk] = engine->mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1];
        }
        for (;kk<N-1;kk++) {
            y = (engine->mt[kk]&UPPER_MASK)|(engine->mt[kk+1]&LOWER_MASK);
            engine->mt[kk] = engine->mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1];
        }
        y = (engine->mt[N-1]&UPPER_MASK)|(engine->mt[0]&LOWER_MASK);
        engine->mt[N-1] = engine->mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1];
        engine->mti = 0;
    }
    y = engine->mt[engine->mti++];
    /* Tempering */
    y ^= (y >> 11);
    y ^= (y << 7) & 0x9d2c5680;
    y ^= (y << 15) & 0xefc60000;
    y ^= (y >> 18);
    return y;
}

unsigned mt19937_randuint32(mt19937_engine *engine){
    return mt19937_randuint32_nolock(engine);;
}

/* generates a random number on [0,0x7fffffff]-interval */
int mt19937_randint31(mt19937_engine *engine)
{
    return (int)(mt19937_randuint32(engine)>>1);
}

/* generates a random number on [0,1]-real-interval */
double mt19937_randreal1(mt19937_engine *engine)
{
    return mt19937_randuint32(engine)*(1.0/4294967295.0);
    /* divided by 2^32-1 */
}

/* generates a random number on [0,1)-real-interval */
double mt19937_randreal2(mt19937_engine *engine)
{
    return mt19937_randuint32(engine)*(1.0/4294967296.0);
    /* divided by 2^32 */
}

/* generates a random number on (0,1)-real-interval */
double mt19937_randreal3(mt19937_engine *engine)
{
    return (((double)mt19937_randuint32(engine)) + 0.5)*(1.0/4294967296.0);
    /* divided by 2^32 */
}

/* generates a random number on [0,1) with 53-bit resolution*/
double mt19937_randres53(mt19937_engine *engine)
{
    unsigned a,b;
    a=mt19937_randuint32_nolock(engine)>>5;
    b=mt19937_randuint32_nolock(engine)>>6;
    return(a*67108864.0+b)*(1.0/9007199254740992.0);
}
/* These real versions are due to Isaku Wada, 2002/01/09 added */

/* initializes mt[NN] with a seed */
mt19937_64_engine *create_mt19937_64(unsigned long long seed)
{
    unsigned mti;
    mt19937_64_engine *engine=(mt19937_64_engine*)malloc(sizeof(mt19937_64_engine));
    engine->mt[0] = seed;
    for (mti=1; mti<NN; mti++){
        engine->mt[mti] =  (6364136223846793005ULL * (engine->mt[mti-1] ^ (engine->mt[mti-1] >> 62)) + mti);
    }
    engine->mti=mti;
    return engine;
}

/* initialize by an array with array-length */
/* init_key is the array for initializing keys */
/* key_length is its length */
mt19937_64_engine *create_mt19937_64_by_array(const unsigned long long *init_key,
             size_t key_length)
{
    unsigned long long i=1, j=0, k= (NN>key_length ? NN : key_length);
    mt19937_64_engine *engine=create_mt19937_64(19650218ULL);
    for (; k; k--) {
        engine->mt[i] = (engine->mt[i] ^ ((engine->mt[i-1] ^ (engine->mt[i-1] >> 62)) * 3935559000370003845ULL))
          + init_key[j] + j; /* non linear */
        i++; j++;
        if (i>=NN) { engine->mt[0] = engine->mt[NN-1]; i=1; }
        if (j>=key_length) j=0;
    }
    for (k=NN-1; k; k--) {
        engine->mt[i] = (engine->mt[i] ^ ((engine->mt[i-1] ^ (engine->mt[i-1] >> 62)) * 2862933555777941757ULL))
          - i; /* non linear */
        i++;
        if (i>=NN) { engine->mt[0] = engine->mt[NN-1]; i=1; }
    }

    engine->mt[0] = 1ULL << 63; /* MSB is 1; assuring non-zero initial array */
    return engine;
}

/* generates a random number on [0, 2^64-1]-interval */
static unsigned long long mt19937_64_randuint64_nolock(mt19937_64_engine *engine)
{
    unsigned long long i;
    unsigned long long x;
    static const unsigned long long mag01[2]={0ULL, MATRIX_A_64};

    if (engine->mti >= NN) { /* generate NN words at one time */

        /* if init_genrand64() has not been called, */
        /* a default initial seed is used     */
        for (i=0;i<NN-MM;i++) {
            x = (engine->mt[i]&UM)|(engine->mt[i+1]&LM);
            engine->mt[i] = engine->mt[i+MM] ^ (x>>1) ^ mag01[x&1ULL];
        }
        for (;i<NN-1;i++) {
            x = (engine->mt[i]&UM)|(engine->mt[i+1]&LM);
            engine->mt[i] = engine->mt[i+(MM-NN)] ^ (x>>1) ^ mag01[x&1ULL];
        }
        x = (engine->mt[NN-1]&UM)|(engine->mt[0]&LM);
        engine->mt[NN-1] = engine->mt[MM-1] ^ (x>>1) ^ mag01[x&1ULL];

        engine->mti = 0;
    }

    x = engine->mt[engine->mti++];

    x ^= (x >> 29) & 0x5555555555555555ULL;
    x ^= (x << 17) & 0x71D67FFFEDA60000ULL;
    x ^= (x << 37) & 0xFFF7EEE000000000ULL;
    x ^= (x >> 43);
    return x;
}

unsigned long long mt19937_64_randuint64(mt19937_64_engine *engine){
    return mt19937_64_randuint64_nolock(engine);
}

/* generates a random number on [0, 2^63-1]-interval */
long long mt19937_64_randint63(mt19937_64_engine *engine)
{
    return (long long)(mt19937_64_randuint64(engine) >> 1);
}

/* generates a random number on [0,1]-real-interval */
double mt19937_64_randreal1(mt19937_64_engine *engine)
{
    return (mt19937_64_randuint64(engine) >> 11) * (1.0/9007199254740991.0);
}

/* generates a random number on [0,1)-real-interval */
double mt19937_64_randreal2(mt19937_64_engine *engine)
{
    return (mt19937_64_randuint64(engine) >> 11) * (1.0/9007199254740992.0);
}

/* generates a random number on (0,1)-real-interval */
double mt19937_64_randreal3(mt19937_64_engine *engine)
{
    return ((mt19937_64_randuint64(engine) >> 12) + 0.5) * (1.0/4503599627370496.0);
}

void destroy_mt19937_64(mt19937_64_engine *engine){
    free(engine);
}
